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The process of protein synthesis in biological systems resembles a one dimensional driven lattice 
gas in which the particles have spatial extent, covering more than one lattice site. We expand the 
well studied Totally Asymmetric Exclusion Process (TASEP), in which particles typically cover a 
single lattice site, to include cases with extended objects. Exact solutions can be determined for 
a uniform closed system. We analyze the uniform open system through two approaches. First, a 
continuum limit produces a modified diffusion equation for particle density profiles. Second, an 
extremal principle based on domain wall theory accurately predicts the phase diagram and currents 
in each phase. Finally, we briefly consider approximate approaches to a non-uniform open system 
with quenched disorder in the particle hopping rates and compare these approaches with Monte 
Carlo simulations. 
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Introduction 



The process of protein synthesis, important in biolog- 
ical systems, has been the focus of intensive study over 
the last few decades. We concentrate on protein synthe- 
sis in prokaryotes, particularly Escherichia coli, because 
it is relatively simple and well studied. The mechanism 
consists of ribosomes "reading" the codons of messenger 
RNA (mRNA) as the ribosomes move along an mRNA 
chain, and the recruitment and assembly of amino acids 
(appropriate to the codons being read) to form a pro- 
tein. (See, e.g., Q, for more details.) Known as "trans- 
lation," this process is often described as three steps: 
initiation, where ribosomes attach themselves, one at a 
time, at the "start" end of the mRNA; elongation, where 
the ribosomes move down the chain in a series of steps; 
and termination, where they detach at the "stop" codon. 
Since ribosomes cannot overlap, their dynamics is subject 
to the "excluded volume constraint." Apart from being 
impeded by another ribosome (steric hinderance), a ri- 
bosome cannot move until the arrival of an appropriate 
transfer RNA, carrying the appropriate amino acid (a 
combination known as aminoacyl-tRNA, or aa-tRNA). 
Thus, the relative abundances of the approximately 60 
types |2] of aa-tRNA have significant effects on the elon- 
gation rate. Assuming reactant availabilities in a cell are 
in their steady state, with a time-independent concentra- 
tion of ribosomes and aa-tRNA, there would be an ap- 
proximately steady (average) current of ribosomes mov- 
ing along the mRNA, resulting in a specific production 
rate of this particular protein. Our goal is the predic- 
tion of the protein production rates for various mRNA's, 
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as a function of the concentration of ribosomes and aa- 
tRNA's. 

The process of translation is well suited to modeling 
using a driven lattice gas in one dimension. In most 
relevant studies of one dimensional driven lattice gases, 
particles are injected at some rate on one end of a chain 
of discrete lattice sites, then hop down the chain one 
site at a time with another rate, and finally exit the 
chain at the other end with a third rate. These three 
rates correspond to the rates of initiation, elongation, 
and termination. The excluded volume constraint is im- 
plemented by insuring that each site can accommodate 
at most one particle. Because the dynamics is stochastic, 
even this "simple" model, though soluable exactly |^ , 
is already not trivial. However, we believe two essen- 
tial aspects of translation are missing from this model. 
First, if we model a codon by a lattice site, the ribo- 
some would cover, typically, a dozen sites 0,0 ■ Second, 
there is non-uniformity in the hopping (elongation) rates 
along the chain, because a ribosome has to "wait" for 
the appropriate aa-tRNA before continuing, and the rel- 
ative abundance of the different aa-tRNA's is far from 
unity. Remarkably, the first issue was explored as early 
as 1968 0,0 J though only at the deterministic, "mean 
field" level. 

In this paper, we present studies which address both of 
these issues, extending the work on the "simple model" 
known as TASEP (namely, single-site coverage. To- 
tally Asymmetric Simple Exclusion Processes with open 
boundaries) 0. Our methods involve both Monte Carlo 
simulations and modern analysis techniques, including 
domain wall theory [Tol |. We have confirmed all key re- 
sults from earlier studies 0, , and several new insights 
have emerged. The paper is organized as follows. The 
details of the model are delineated first, along with brief 
summaries of known results. Section ^1 is devoted to a 
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dosed (i.e., periodic) system with particles of arbitrary 
size. Though not a direct model for translation, TASEP 
on a "ring" provides simple solutions as well as useful 
insights in the form of exact relations for relevant pa- 
rameters, such as current-density relations. In Section 
mil we turn to the central topic: TASEP with extended 
objects and open boundaries. Section IIVI is devoted to 
non-uniform hopping (elongation) rates. We close with 
a brief summary and speculate on the relevance of our 
model as a mechanism for the nonlinear relationship be- 
tween mRNA ar id p rotein levels observed in biological 
experiments [Til Il2j| . 

I. THE DRIVEN LATTICE GAS AS A MODEL 
FOR PROTEIN SYNTHESIS 

A. Model specifications 

We model an mRNA with N codons as a chain of sites, 
each of which is labeled by i. The first and last sites, 
i = 1,N, are associated with the start and stop codons 
respectively. At any time, attached to the mRNA are 
M ribosomes (also referred to as "particles"), which we 
label by a (a = 1,...,M). Being a large complex of 
molecules, each ribosome will cover i sites (codons), with 
i = 12 typically By contrast, nearly all studies of 

the asymmetric simple exclusion processes (ASEP) are 
devoted to ^ = 1. Any site may be covered by a single 
ribosome or none. In case of the latter, we will refer to the 
site as "empty" or "occupied by a hole." For convenience, 
we define M as the number of holes on the chain, so that 

M + iM^N. (1) 

For open systems, a ribosome at the end can be attached 
without covering all i codons, so that this equality is only 
approximately true. To locate the ribosome, we arbitrar- 
ily choose the lowest site covered. For example, if the 
first i sites are empty, a ribosome can bind in an ini- 
tiation step, and then it is said to be "on site i — 1." 
Therefore, a complete specification of the configuration 
(state, or micro-state) of the mRNA is the set of loca- 
tions: {ia}- One disadvantage of this labeling is that, 
with each initiation event, the a of every ribosome will 
change (increase by unity) . Alternatively, we can use site 
occupation numbers 

J 1 if site i is covered by any part of a ribosome 
' ~ |_ if site i is empty . 

With these conventions, we define several density pa- 
rameters: 

• Pr = M/N is the ribosome (or particle) density; 

• p = M£/N = J2i "n-i/N is the coverage density; 

• p/i = 1 — p is the hole density; and 



• ps = Pr + Ph is defined for convenience. 



All of these quantities are time dependent, because ini- 
tiation and termination occur independently. For math- 
ematical reasons, we will first consider a closed system 
(with periodic boundary conditions, i.e., the ends of the 
chain tied to form a ring), for which these densities are 
fixed and Eqn. Q holds strictly. As will be clear later, 
it is also convenient to label configurations by specifying 
the number of holes between successive ribosomes: {ha}, 
where ha is the number of holes in front of the a*'' ri- 
bosome. In the terminology of traffic models, ha is also 
known as the "headway" of this particle. Though not 
absolutely necessary, we could define ho as the number 
of holes behind the first ribosome. 

Next, we specify the dynamics of our model. An at- 
tached ribosome located at site i will move to the next 
site (i -|- 1) with a rate ki, provided site i -\- £ is empty. 
For Monte Carlo simulations, it is convenient to update 
configurations in discrete time units. Then, it is better to 
use probabilities Pi (0 < pi < 1), so that a ribosome on 
site i will be moved or not with probability pi or I — pi, 
respectively. We purposefully associate these hopping 
probabilities with a site because a site is associated with 
a particular codon. Thus, the jump rate from that site 
may depend on the relative abundance of the appropriate 
aa-tRNA. Apart from these probabilities, another aspect 
of our stochastics is random sequential updating: i.e., 
during each Monte Carlo step (MCS), M + 1 particles 
are chosen at random, in sequence, to attempt moves. 
They are selected from a pool that includes the M par- 
ticles on the lattice plus another unbound particle that 
can initiate if there are £ holes at the beginning of the 
chain. Let us illustrate with a few examples. First, po is 
associated with the start codon and, if the first £ sites 
are empty, a particle will be placed on the i = 1 site with 
this probability. Next, a random particle (say at site i) 
is chosen and, provided it has a hole in front (n^+f = 0), 
will be moved with probability pi . Naturally, it will not 
be necessary to check for headway for the "last" ribosome 
{a = M ). Finally, the stop codon will be associated with 
Pn. For simulations of the closed system, there will be no 
"beginning" or "end," so that there are no special steps 
for initiation or termination. 

In our computational studies, 100 identical systems of 
N sites are simulated in parallel to obtain good statistics. 
Simulations of closed systems begin with particles evenly 
distributed around the ring and run for 3600 MCS to 
ensure that steady state is reached. Open systems begin 
empty and are run for 12,000 MCS (for N < 500) or lOOA^ 
MCS (for N > 500) to reach steady state. After steady 
state is attained, data including the current and density 
distribution can be collected. Density data is typically 
collected everylOO MCS. We often use continuous time 
Monte Carlo [I3 because it runs far more quickly than 
and provides the same results as standard Monte Carlo. 
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B. Brief survey of known results 

Extensive investigations of the simple totally asymmet- 
ric exclusion process (TASEP, defined as point particles 
hopping with unit rate along a line) with open bound- 
aries can be found in the literature. Simulations have 
been performed [Ijl, and exact analytic results for the 
steady state exist |3, Q . Depending on the initiation (or 
injection) and termination (or depletion) rates, the sys- 
tem will settle into one of three phases. Introduced above 
as po and pN respectively, the initiation and termination 
probabilities are mostly referred to as simply a and /3 in 
the literature. From their dominant characteristics, the 
three phases are known as low density, high density, and 
maximal current. A phase diagram in this a-P plane has 
been determined, showing second order transitions be- 
tween the maximal current phase and the others, as well 
as a first order transition between the high- and low- 
density regions. Subtle correlations further divide the 
last two into subregions. When disorder is introduced, 
i.e., not all the p^'s are equal, then methods for exact 
analytic approaches fail (except in the extremely dilute 
limit, where only the motion of a single particle is of 
concern p^l. Indeed, even a single slow rate in a closed 
system poses serious difficulties |16i, 13 • However, 
Kolomeisky obtained approximate steady state so- 
lutions and phase diagrams for an open system with a 
single nonuniform rate in the bulk by splitting the sys- 
tem into two smaller systems connected by the nonuni- 
form rate. Tripathy and Barma considered a closed 
system, but with a finite fraction of identical slow sites. 
Based on a combination of Monte Carlo simulations and 
numerical solutions of mean field equations, they found 
current-density relations. Though not a model for elon- 
gation, a related problem is "particlewise disorder," in 
which the unequal hopping rates are associated with the 
particles rather than the sites 21\. Further references 
on TASEP with disorder may be found in a recent re- 
view |2.2J . Also indirectly related to our one dimensional 
models are driven lattice gases with quenched disorder 
in higher dimensions . Finally, we mention that 

there are many studies on the ASEP in which a particle 
has a finite probability of stepping backwards z5\ . Back 
steps are not generally believed to occur in elongation, 
and we will not consider such processes. All of these 
studies are restricted to £ = 1. 

Systems with extended objects (^ > 1) have been 
rarely investigated, despite their introduction over three 
decades ago as a model for biopolymerization . Using 
a mean field approach, MacDonald et al. set up mean 
field equations for the average site occupation {rii) and 
considered both closed W\ and open S| systems. In the 
former case, exact solutions were found, leading them to 
a current vs. density relation. For the latter, the authors 
resorted to numerical solutions to find the phase diagram 
for a variety of initiation and termination rates. A phase 
diagram similar to the simple TASEP, as well as non- 
trivial density profiles and the associated currents, was 



obtained. More recently, there is renewed interest in this 
problem. Naming this system TASEP," Sasamoto and 
Wadati focused on the time dependence of M par- 
ticles in an infinite lattice and, using the Bethe ansatz, 
found exact results for the conditional probability that 
the particles are found at certain sites given an initial 
configuration. The main conclusion is that the dynamics 
of TASEP lies in the same universality class as ordi- 
nary TASEP. This line of inquiry has been further gen- 
eralized to a system containing a distribution of particle 
sizes [23, 13 • Though these investigations produced in- 
teresting results, they are not applicable to our situation, 
namely, finite systems with open boundaries. In particu- 
lar, for finite lattices, these studies are restricted to closed 
systems, for which the stationary states are trivial, as we 
will recapitulate in Section^ Finally, a recent work by 
Lakatos and Chou considered uniform open systems 
with extended objects. Using a discrete Tonks gas parti- 
tion function, they derived the current vs. density rela- 
tion first presented by MacDonald et al. . Via a refined 
mean field theory, they extended this result to predict 
currents and bulk densities for the open system, which 
they confirmed by Monte Carlo simulations. The phase 
diagram and its properties are consistent with those ini- 
tially obtained by MacDonald and Gibbs Q . We are not 
aware of published results on open systems with both 
extended objects and quenched disorder. 



II. TASEP OF EXTENDED OBJECTS ON A 
RING 

For simplicity and mathematical reasons, it is conve- 
nient to discuss a uniform closed system with periodic 
boundary conditions. Here, M identical particles of size 
£ move on a ring of N sites. Due to the excluded vol- 
ume constraint, there are M ^ N — IM uncovered sites 
(holes). The system evolves by discrete time steps, with 
random sequential updates. Specifically, during each 
MCS, M particles are chosen at random in sequence, and 
each is moved forward one site provided there is head- 
way. Since we have a stochastic dynamics, the complete 
description of this system is P (C, t) , the probability that 
it is found in configuration C after t steps (starting from 
some initial Co). To label a configuration, we choose an 
arbitrary particle to be the first (a = 1) and supply {/la}, 
the set of number of holes in front of the a*'' particle 
(with a = 1, M). Clearly, we may also think of a con- 
figuration as a series of M "gaps" (between the particles) 
with ha being the number of holes in the a*'* gap. So, 
P {hi, h2, ■■■hM',t) is an explicit form for P{C,t). Note 
that, since the system is closed, there is a constraint on 
ha, i.e., 

a. 

is a constant in time. 
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Random sequential updating can be translated into an 
equation which governs the time evolution of P (C, t) , 
namely a master equation. Starting with the ini- 
tial P(C, 0) — 6{C,Ca) (where 6 is the Kronecker 
delta), P{C,t) is expected to settle into a unique, time- 
independent distribution, P* (C), which we will refer to 
as the "steady state." If this system were evolving to- 
wards thermal e quil ibrium, the dynamics would satisfy 
detailed balance j^, and P* (C) would be given by the 
well-known Boltzmann factor. However, the dynamics of 
our system definitely violates detailed balance, so that, 
associated with a non- equilibrium steady state, P* (C) 
is not known in general. Fortunately, this closed system 
belongs to a class for which a simple solution is known 
|3lj | ; namely, every configuration occurs with equal prob- 
ability [s^. Thus, P* (C) is precisely the reciprocal of 
the total number of configurations consistent with the 
given parameters {N, M, £) . With such a simple distri- 
bution, we can compute many quantities of interest, such 
as the probability distribution of the current (and hence 
the average current) and the headway. 

Apart from overall factors, the total number of config- 
urations is just Z{M, M), the total number of lists {ha} 
subject to Eqn. Because this is a well-known com- 
binatorial problem (appearing in, e.g., the Bose gas), we 
simply quote the result: 



Z{M,M) = 



{M + M- 1)! _ /M + M - 1 
M!(M- 1)! " V M 



The actual number of distinct configurations is 
(N/M) Z{M,M), because there are N lattice sites on 
which the first particle can be placed but the AI parti- 
cles are identical. Thus, 



P* (C) 



MlMl 



N{M + M - 1)1 



independent of C. 

To compute the probability distribution of the current 
and its average, we need more detailed information on 
the above "partition." Defining the current, J, as the 
number of particles which moved in one step (normalized 
by the system size N), we see that we will need H, the 
number of gaps with one or more holes. The reason is 
that, for each such gap, the ribosome behind it can move 
and contributes one "unit" to the current, so that 



Since all configurations are equally probable, the statisti- 
cal weight associated with this J is just the total number 
of configurations with a given H. This quantity, denoted 
by Z{H; M, M) in analog to Z(M, Af ), may be found 
from its definition: 



Z{H] M, M) 




where the sum is over all possible lists {ha] and the Kro- 
necker (J's select only those which satisfy Eqn. (JJ and 
have H gaps with ha > 0. Once Z{H; M, M) is known, 
the full distribution for the current is 



p{J-M,M) 



Z{H;M,M) 
Z{M,M) 



Now, the explicit form of Z{H; M, M)) can be obtained 
either through the generating function 



M,H 



1 



1-C 



M 



or by standard combinatorial techniques: 

Thus, the explicit current distribution is 

'M\ /M - 1\ / /M + M -1 



P{J) = 



H \H-l 




M 



(2) 



where H on the right stands for JN. An alternate form, 
showing the dependence of this distribution on the con- 
trol parameters (M, M,N = M + £M) is 



p( J I M, M, TV) 



1 



MM M + M 



JN 



(M- JiV)!(A'f- JiV)! 
MlMl''^ 



{JNy. 



To illustrate this distribution, we show in Figure Q both 
this prediction and simulation data for the case of TV = 
200, M = 15, and £ = 12. Clearly, there is excellent 
agreement between theory and simulation. 

Containing less information, but easier to grasp, is the 
average current J = ^ Jp (J). Its computation is some- 
what easier, since it is Y.h=i HZ{H; M, M)/NZ{M.M), 
with the numerator easily gleaned from dnWuiC, ■ 
The result is 



M 



M 



N M + M -1 



(3) 
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FIG. 1: Distribution of currents p{J) for N = 200, M = 15, 
and I = 12. Squares (connected by dashed lines) are val- 
ues predicted by Eqn. and triangles are values observed 
in Monte Carlo simulations. A single lattice was simulated 
to steady state, and instantaneous current (H/N) was deter- 
mined every 100 MCS thereafter for 1.2 x 10*^ MCS. 



As we will see, the dependence of this average current 
on the density of particles plays a central role. Ex- 
pressing this quantity in terms of pr = M/N, we have 
J = (Or (1 — ^Pr) / [1 ^ ~ 1) — 1/^] ■ An appcahug 
form, which displays both the intensive nature of J and 
its underlying particle-hole symmetry, is 



J 



PrPh 



Ps - l/N 



(4) 



A third form, frequently referred to in the literature as 
the "current-density relationship," is writing J as a func- 
tion of p, the coverage density [p e [0, 1]): 



J{P) 



1 



II- p + p/l-l/N 



(5) 



In the limit iV ^ oo, this result was first presented by 
MacDonald et al. (3 • A generalization of the well-known 
expression for ^ = 1 (i.e., J = p{l — p)), this J{p) is no 
longer symmetric about p = 1/2. Instead, the optimal 

density increases from 1/2 to VI y/ ^1 + , while the 

maximum current is lowered from 1/4 to ^1 + 

As these quantities will appear frequently, we will denote 
them by 



i + VI 



and J ; 



i + Vi 



(6) 



To appreciate this shift graphically, we present, in Figure 
El both analytic and simulation results. In addition, to 
show the effects of the finite size corrections (due to the 
1/A^ term), we include a curve of the limiting form |^|29|| 



J 



PrPh 



(7) 
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FIG. 2: Current J vs. coverage density p — M£/N for closed 
systems. Symbols are Monte Carlo results, solid curves are 
predicted values from Eqn. Q, and broken curve is prediction 
from Eqn. 0. Triangles are for £ = 2, A'' = 40, x 's for £ = 5, 
N = 80, and squares for £ = 12, N = 150. J was determined 
by averaging over 1.2 x 10^ MCS and 100 identical systems 
after steady state was reached. 



for the iV = 40 case. 

In connection with these expressions, a conclusion for 
conditional probabilities can be drawn. Since J is pre- 
cisely the joint probability of finding a "covered" -hole 
pair, we see that Pr/ps is the probability that site i is cov- 
ered, given that site i-|-l is empty and, similarly, Phl Ps is 
the probability that site i -I- 1 is empty, given that site i is 
covered. These conditional probabilities will play a role 
in our understanding of the behavior of open systems. 

In addition to the average current, we can also compute 
its fluctuations exactly. 



Aj2 = ^j2p(j)_ J2 



PsN 



Typical of non-critical thermodynamic systems, in which 
the (fractional) deviations are 0(Af^^/^), the full dis- 
tribution p{J) approaches the standard Gaussian form: 



exp 



fp. {J/ J -I) 



Indeed, this is the form we see 
in the example shown in Figure ^a-bove. 

Finally, we turn to another quantity of interest: the 
statistics of "headway." Since there are M gaps in each 
configuration and there are Z{M, M) configurations, we 
have a total of MZ{M, M) gaps. Out of these, we wish 
to compute the number of gaps which contain precisely 
TO holes, which we denote by Z(to; M, M). From its def- 
inition 



Z{m- M, M)= ^ (5 Af , ^ /i„ (to, K) 

{ha} \ a / \ a 
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FIG. 3: Distribution of head spacings p(m) for — 100, 
M — 6, and £ — 12. Curve is prediction from Eqn. and 
squares are simulation values. 100 lattices were simulated to 
steady state, and head spacing m for any particles on the 
first lattice sites was determined every 100 MCS thereafter 
for 1.2 X 10^ MCS. 



we find the associated generating function 



M,m 



= M 



1 



1 



M-1 



leading to the distribution of head spacings 



p{m; M, M) 



Z{m;M,M) 
MZ{M, M) 

fM+M-m-2\ 
\ M-m I 



( 



M+M- 
M 



(8) 



In the hmit of large iV, this expression simplifies to 
p{m) — — 

Ps \Ps 

This distribution is reproduced faithfully in Monte Carlo 
simulations, as shown by an example in Figure |3| It 
is easy to understand this result intuitively if we re- 
gard Ph/ Ps as the probability of having a single hole in 
the headway. With independent hole statistics, we have 
p{m) cx (ph/ps) ■ The average number of holes in the 
headway and the associated standard deviation can be 

computed easily: ph / Pr and ^J {ph/ Prf' + Ph/ Pr, respec- 
tively. 

Though the system considered in this section, parti- 
cles traveling on a ring, bears little resemblance to the 
translation process, it is sufficiently simple for us to de- 
rive a number of exact results. Apart from their own in- 
terest, these results provide crucial insights, such as the 



current-density relationship, for postulating appropriate 
equations in a coarse grained, mean field approach to the 
physical problem at hand. 



III. EXTENDED OBJECTS IN OPEN SYSTEMS 

In this section, we attempt the next step towards a 
realistic model for protein synthesis by considering sys- 
tems with open boundaries. The section is organized as 
follows. We describe the problem and introduce some 
terminology. Next, we study the effects of the open 
boundaries alone, by analyzing the continuum limit of 
a symmetric exclusion process. Finally, we consider the 
asymmetric exclusion process with open boundaries. We 
present the phase diagram found from simulations and 
an extremal principle analysis, and we show the ability 
of the continuum limit to predict steady state density 
profiles. 

In the open system, the first site {i ~ 1) is no longer 
"in front of" the i = N site. Instead, a particle (of extent 
£) will be placed at i = 1 with probability a (previously 
labeled by po), provided all the first £ sites are empty. 
This models the initiation process. As for elongation, we 
continue to restrict ourselves to uniform rates here. In 
the language of TASEP, every (randomly) chosen particle 
will move by one site (increasing i by unity) if it has some 
headway. For particles on the last £ sites, there is no 
hinderance, so they will always move if chosen. Finally, 
to simulate termination, a particle on the N^^ site will be 
removed from the system with rate f3 (previously labeled 
by pn). To repeat, in a Monte Carlo step (MCS), M + 1 
particles are chosen randomly (in sequence) to attempt 
a move. They are selected from a pool including the M 
particles on the lattice plus an unbound particle that may 
initiate. 

Since this system is no longer closed, M and M are 
fluctuating quantities. Of course, their average values 
will be controlled by the rates a and (3. Our goal is to 
find the average densities and the average current of such 
a model, as functions of (a,/3). In the £ = 1 case, it is 
known that the system exhibits three different "phases" 
as a and /3 are varied, only one of which resembles the 
closed system above (in the sense that the current ap- 
proaches J for large N). Beyond overall averages, we seek 
the density profile, which will not only be non-trivial, 
but which also displays drastically different properties 
as we move about in the a-/? plane. The main goal of 
this paper is to study the effects of ^ > 1 on both the 
phase diagram and these density profiles. Unfortunately, 
P* (C) for a system with open boundaries is not known 
in general. Even for the £ = 1 case, only a limited set of 
quantities may be computed exactly. To make progress, 
we resort here to a more phenomenological approach, in 
the spirit of hydrodynamics or Lanudau-Ginzburg free- 
energy functionals. Considering coarse-grained densities 
and the continuum limit, we postulate equations of mo- 
tion, based on some of the properties of the closed sys- 
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tern. In principle, such equations can be "derived" from 
the master equation for P (C, t) , using the mean- field ap- 
proximation (i.e., ignoring all correlations). In a sense, 
this is also the approach of MacDonald et al. i2| , except 

that they focused on n-^^ (t) and nl"^ (i) , the probabil- 
ities for site i to be occupied by a ribosome and a hole, 
respectively, at time t. By keeping the spatial co-ordinate 
discrete, they faced difference equations and succeeded in 
finding solutions only numerically. In contrast, by using a 
continuous spatial co-ordinate, we have ordinary differen- 
tial equations instead, giving us analytic solutions. How- 
ever, lost in the notion of "coarse-graining" are the period 
£ structures which feature prominently behind "block- 
ages." Nevertheless, our approach appears to capture 
the essential (gross) features of these systems. 

In our heuristic approach, we imagine N to be large 
enough to justify taking the continuum limit, i.e., replac- 
ing the discrete site label z by a continuous co-ordinate: 
X. For simplicity, define 

X = i/N 

so that a; lies within the unit interval. (If physical units 
of length are desired, we may introduce a as the lat- 
tice spacing, corresponding to the length of a codon, i.e., 
three bases of the mRNA. Then L^rna = Na would 
be the length of the mRNA in question.) Similarly, con- 
tinuous local densities will take the place of the discrete 
occupation variables. For example, the coverage density 
p{x) will be used instead of Ui. Since the maximum oc- 
cupancy is unity, the hole density is just 

Ph{x) = l- p (x) . 

Following our considerations above, we also define the 
ribosome density by 

Pr (x) = p{x) II . 

Note that, despite the continuum limit, these equations 
display the meaning of ^, which serves as a measure of 
the "size" (or "extent" ) of ribosome. Of course, we are 
also interested in their time dependence, so that we must 
consider p {x, t) in general. Now, ribosomes rarely de- 
tach from the mRNA during elongation, so that we are 
justified in regarding these densities as conserved fields. 
Thus, the appropriate equation of motion is the continu- 
ity equation, i.e., 



-V • Jr 



(9) 



where Jr {x, t) is the local (ribosome) current. Our first 
task is to find how this current depends on the (lo- 
cal) density pr {x^t), i.e., to find the functional form for 
Jr [pr] ■ Then we will arrive at an acceptable equation for 
Pr- To find the steady state profile p* {x), which satis- 
fies dtPr = 0, we see that this state is associated with 
a constant (x, t-independent) current J*. Therefore, our 
problem consists of finding the solution to 



subject to the appropriate boundary conditions. To keep 
the notation simple, from here on we will drop the sub- 
script r and write J or J {x,t) for the local current (as 
well as J* for the constant current in the stationary 
state). 



A. The case of free diffusion and an effective 
density variable for extended objects 

To show how we build an appropriate equation for the 
TASEP, let us begin with a study of the effects of open 
boundaries alone. In other words, let us consider a sym- 
metric exclusion process for extended objects, i.e., a sys- 
tem of large particles diffusing freely on a line, subjected 
to the excluded volume constraint and the controls a, (3. 
For simulations, a randomly chosen particle is moved one 
site forward or backward with equal probability (0.5), 
provided it does not run into its neighbor. The only 
exception is the first particle, which is prohibited from 
jumping backwards into the "source." 

For the i — 1 case {p = pr), the steady state profile is 
trivially linear: 

p* {x) = {l^x)p* {0) + {x)p*{l) 

with current 

J*-^[P*(0)-P*(1)] . 

(Note that, since the current is controlled by the gradient 
of the local density only, it vanishes necessarily in the 
N ^ oo limit. So, we must keep the N explicitly here.) 
Meanwhile, the boundary densities are fixed by matching 
the injection/depletion rates to the internal current: 

apliO) = a[l - p* m = J* ^ f3p* {!) , 

and our problem is completely solved. These well known 
results can be traced to the fact that Eqn. assumes 
the form of the simple diffusion equation dtp oa d^p, since 
J [p] (X — Vp. Though there are many ways to arrive at 
this result, it is less clear how to generalize it to the case 
of extended objects. In particular, as displayed in Figure 
21 the profile from simulations (for N = 200, ^ = 12,a = 
/9 = 1) is far from linear! Here, we will show that there 
is a natural generalization for the (particle) density in 
the case of extended objects and that its profile is again 
linear. 

We proceed by returning to the basics, starting with 
the current Jr being proportional to both the conductiv- 
ity (mobility) and the drive. For the former, we take the 
result from the previous section. There, the drive is con- 
stant, so that the right hand side of Eqn. (0, prPh/ps, 
can be interpreted as the conductivity or mobility (ne- 
glecting the finite-size effect term from the closed sys- 
tem). Meanwhile, the drive for free diffusion should be 
the gradient of a pressure V, so that 



Jr [Pr 



constant 



Jr [p] = D 



PrPh 



Ps 



(10) 
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where D is a constant, to be fitted to data. Since we 
have scaled the system size to unity, we should keep in 
mind that _D is a quantity of 0(l/iV), to be consistent 
with the continuum limit. For the pressure V ^ we follow 
the standard route of statistical mechanics and write 



V 



5H 

Spr 



where 7i is a free energy functional (e.g., the Landau- 
Ginzburg "Hamiltonian" in case of the ordinary lattice 
gas). Here, we have a non-interacting system, so that a 
reasonable form for H is just the entropy ,35, ] : 



dx [pr In pr + ph In Ph - Ps In Ps 



Using 



p/i = 1 - £pr 



Ps 



l-i£~l)pr 



and carrying out the steps, we arrive at a modified diffu- 
sion equation: 



dtPr = Ddx [Ps ^d^pr] 



(11) 



Note that, for £ > 1, the effective "diffusion constant," 
D/ p^, is density- dependent. Thus, the steady state profile 
will not be linear in x. Instead, it satisfies 



D dp; 
(Ptf 9x 



= -J* 



(12) 



which will definitely lead to non-linear profiles. 

Although we have derived this diffusion equation via 
phenomenological techniques, wc speculate that it could 
be obtained through a more rigorous derivation. Specifi- 
cally, if a matrix product description of the £ > 1 sy stem 
becomes available, the methods of Derrida et al. |3p| may 
extend to calculate a large deviation functional, which 
would play the role of a nonequilibrium free energy for 
this system. Should such a functional exist, we expect 
that our diffusion equation would result from a variation 
of this functional. 

To find p* explicitly, notice that, due to = 1 — 
{£ — 1) Pr, the left hand side is just a simple derivative of 
Thus, the profile (x) will be linear. However, 
Ps does not reduce to a sensible variable as ^ ^ 1, leading 
us to define 



1 



£-{£--l)p* 



Not only does x li^ in the unit interval; it reduces to 
the usual particle density p* in the limit £ ^ 1. Most 
importantly, it is uniquely related to p* for all £. First 
introduced in (2^, x is the "natural" generalization of 
the particle density for extended objects. We will refer 
to it as the effective particle density (EPD). Meanwhile, 
a hole remains a single-site entity, so its density needs no 



modification. In terms of the EPD, the various densities 
are 



Ps = 

Pr = 
Ph ^ 



i + (.£~i)x' 

X 

i + i£-i)x' 
i-x 

l + i£-l)x 



SO that X is just 



X = PriPs- 



(13) 



The current-density relationship (Eqn. [7J for the steady 
state in the ring is again a simple product: 



Jring XPh • 



(14) 



Note that this product will serve as the mobility in our 
mean field approach here and can be expressed in terms 
of X alone: 



XPh 



x(i - x) 

l + i£-l)x 



The main advantage of using the EPD is the re- 
emergence of the familiar combination X (1 ^ x) the 
numerator. Finally, corresponding to the optimal den- 
sity (Eqn.ini), we have 



a quantity which plays a significant role in systems with 
open boundaries. Note that the corresponding hole den- 
sity (1 — /3) also assumes the same value, so that J is just 
X^ . The underlying "particle-hole" symmetry generalizes 
to one under exchange of x <^ , or 



X^ 



(1-x) 



i + (^-i)x 

or, in terms of the more physical p, 

1-p 



p^ 



l-p + p/£ 



We now find that Eqn. (|12() simplifies to the familiar 
equation from ordinary diffusion: 



ox 



-J* 



As a result, the profile in terms of x is again linear in x 
[STf. For completeness, we write the solution: 



with 



X(a;) = (l-x)x(0) + (a;)x(l) , 



x(o)-x(i)- 



9 



^ 




FIG. 4: Density profile for free diffusion in a system with 
iV = 200 and £ = 12 from p = 1 to p = 0. Symbols are 
Monte Carlo data. 100 identical systems were simulated in 
parallel for 10^ MCS to reach steady state, and then den- 
sity profiles were collected every 100 MCS for an additional 
1.2 X lO'' MCS. Curve is predicted density profile. Inset shows 
the effective particle density profile x, constructed from the 
simulation data. The current observed in simulations was 
J* — 0.00252, which is close with the expected value of 
0.0025 = -D [x (0) - X (1)] for D = l/2iV. 



The explicit solution ensues once the constraints of injec- 
tion/depletion rates are imposed. Note that, in general, 
we would have J* oc D, so that the current inherits the 
0{1/N) from D, in contrast to the 0(1) behavior for 
systems with non-zero drive. 

To close, we illustrate how well the theory agrees with 
data by showing the analytic, non-linear profile of p* (x) 
with boundary conditions p(0) = 1, p{l) = 0: 



as a curve in Figure 01 



1-x 



1 - (1 - i/e)x 



B. TASEP in open systems 

We now turn to the other extreme case, where back- 
ward jumps are cornpletely excluded (TASEP), gener- 
alizing the model of Q to extended objects. With open 
boundaries, this model incorporates another essential fea- 
ture of translation. Our interest is again the average cur- 
rent and density profile, as a function of a and /3, the 
rates of initiation and termination. 



1. Phase diagram from extremal principle and simulations 

Although we are unable to generalize the methods of 
0, 3 to arrive at exact solutions for an open system 
with particles covering £ > 1 sites, we find that the 
phase diagram determined by the mean field techniques 
of 0, can be understood also by extending the ex- 
tremal principle, an hypothesis first proposed by Popkov 



and Schiitz 38' and later exploited successfully by oth- 
ers (e.g., ^£,,4^]). In this approach, the open boundaries 
are regarded as connections to reservoirs with appropri- 
ate densities, so that by keeping the same jump rates as 
in the bulk, a and (3 are realized. Defining and p^ as 
the reservoir densities at the initiation and termination 
boundaries, respectively, the extremal principle relates 
the current in the open system to the J[p\ for a closed, 
periodic system with the same bulk dynamics j4l| : 

J _ i max J[p] for p+ < p < p- 
1 min J[p] for p^ < p < pj^- 

Unfortunately, there is no prescription for finding p_ , p^ 
from a, (3 in general. Exploiting results from the exactly 
soluable £ — 1 case (where p- = a and p+=l — (3) and 
from previous studies of the I > 1 case [a , we argue 
in favor of 



(.a 



p-[a) 



l + a{i-l) 



and p+(P) = l-f3. (15) 



These reservoir densities can be understood as follows. 
Recall that, for a closed system, the probability for 
site i — 1 to be filled, given that site i is empty, is 
Pr/Ps — (p/i) / (1 ^ P + p/^)- We now argue that when 
the first site of an open system is empty and will be filled 
with probability a, it can be thought of as being coupled 
to a reservoir of the appropriate density, i.e., p_ such 
that a = {p-/(-) / (1 — P- + P-/^)- Solving for p_ leads 
to the expression above for p_(a). The expression for 
P-i- (/3) is not readily explained by similar arguments, so 
discussion of its origin will be deferred until the follow- 
ing section. An important feature associated with this 
choice is that the current is symmetric under a <^ (3, as, 
observed in simulations. Most importantly, these choices 
lead to phase diagrams in good agreement with Monte 
Carlo data. 

By combining these choices with the extremal princi- 
ple, we find that, although the result is qualitatively sim- 
ilar to the £ = 1 system, there are quantitative changes 
for the (.> \ case. First, there is a shift in the location of 

transition lines, from 1/2 to X = 1/ (l + \/^) (as shown 

in Figure EJ. Then, the current and bulk densities are 
also modified: 



a(l- 



J («,/?) 



I-I-q(^-I) 

l+/3(£-l) 
J 



for a < (3 < X (low density) 
for (3 < a < X (high density) 
for a,(3 >x (max current) , 

(16) 



and 



p(a,/3) 



p_ for a < (3 < X (low density) 

p+ ioY (3 < a < X (high density) (17) 

p for a, (3 > X (max current) . 



As examples of how well these predictions fit simulation 
data, we show plots of J {1,(3) (Figure and p(a,0.1) 
(Figure [TJ. As expected, when a was rate- limiting (low 
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FIG. 5: Phase diagram predicted by the extremal principle. 
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FIG. 6: Dependence of current on /3 for a = 1, ^ = 12, and 
N = 1000. Symbols are simulation results (determined from 
100 systems simulated in parallel for 1.2 x lO** MCS after 
steady state was reached) and curve is the prediction from 
Eqn. (ITCl . 



density phase) , a bulk coverage density of p- {a) was in- 
duced, and when /3 was rate-limiting, a bulk coverage 
density of p+(/3) was induced. Given the success of the 
extremal principle analysis, it would be appropriate to 
search for more fundamental theories from which the ex- 
tremal principle could be derived. 

The transition between low and high density phases 
is clearly first order, displayed as a jump in the bulk 
density plot (Figure [Tj) . As in the I = \ case To^ , do- 
main wall theory can be used to understand our results. 
In particular, for a = /3 < x, we have observed a shock 
front (between the low and high density regions) diffusing 
freely along the lattice. As a result, the average profile 
is linear, as shown in Figure |H1 To appreciate the shock, 
we display in the inset a typical configuration (showing a 
shock) where each ribosome is represented by a dot. In 
addition, we exhibit another aspect of this co-existence 
by sampling the average (coverage) density in the central 



FIG. 7: Dependence of average coverage density p on a for 
l3 = 0.1, I = 12, and iV = 1000. Symbols are simulation 
results (determined from ICQ systems simulated in parallel 
and sampled every 100 MCS for 1.2 x 10* MCS after steady 
state was reached) and curve is the prediction from p_ (a) and 
p+iP) in Eqn. 



10% of a large {N = 1000) lattice. Compiling a histogram 
(closed diamonds in Figure|5J), we find a bimodal distribu- 
tion typical of systems at a first order transition. Peaks 
are expected to correspond to the densities p- {a) and 
p+ (/3) . To connect with local fiuctuations in the "pure" 
systems, we compile similar histograms for a closed sys- 
tem (i.e., 10% of a ring with 1000 sites) with overall den- 
sity set at P-(ck) and p+{l3). For comparison, we show 
a simple linear combination of such distributions in Fig- 
ure|51(open squares). We expect correspondence between 
densities in the window in the open system, which alter- 
nate between approximately p- (a) and p+ (/3), and an 
appropriately weighted average of window densities in the 
two closed systems. The deviations can be understood as 
the contribution of configurations with the shock in the 
sampling window. Roughly, this may occur about 10% 
of the time, which is also the order of magnitude of the 
deviation from the simple average of pure systems. 

While most of the features presented here were known 
to MacDonald et al. QJg and have also been obtained 
by Lakatos and Chou [23, our efforts are to go beyond 
mean field approaches, showing the different perspective 
offered by domain wall theory and the extremal principle. 
Shock effects on the co-existence line are particularly well 
described by this approach. 



2. Steady state profiles from the continuum approach 

Following an understanding of the phase diagram, we 
turn to more details of the system, namely, steady state 
density profiles. Our continuum approach again leads us 
to differential equations. As considerable insight can be 
gained by studying simple differential equations, we be- 
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FIG. 8: Steady state density profile for a = 0.05, f3 = 0.05, 
= 500, £ — 12. Profile was obtained from simulations of 
100 systems run to steady state and then sampled every 100 
MCS for 1.2 X 10* MCS. (Nonlinearities near the termina- 
tion boundary result from ribosomes tending to "pile up" at 
i = N,N - l,N -21, etc. due to /3 < 1 but with decreasing 
correlations as i decreases.) Inset shows a typical configu- 
ration of this system, with the shock front near the center. 
(Each ribosome is represented by a dot.) 




FIG. 9: Probability distribution of densities in the central 
10% of a system with a = 0.1, /? = 0.1, iV = 1000, and 
i — 12. Actual probability distribution (solid diamonds 
with solid curve) was obtained by simulating 100 systems 
to steady state, sampling the central densities every 100 
MCS for 5 X 10* MCS, and compiling a frequency histogram. 
Peaks correspond to the densities p_ (a) and p+ (/3). For 
comparison, a composite density distribution (open squares 
with dashed curve) was assembled from similar frequency his- 
tograms for two closed systems with average densities p- (0.1) 
and p+ (0.1). (Again, N = 1000, £ = 12, and 10% of the sys- 
tem was sampled.) The composite distribution required a lin- 
ear combination of 53% low density and 47% high density to 
minimize the absolute difference from the actual probability 
distribution. See text for details. 



lieve it is worthwhile to devote a section to the continuum 
approach. 

Starting with the case of free diffusion, we extend Eqn. 
(|ll|l to the driven system by simply adding an "Ohmic" 
term to the current: J = aE, where the drive is asso- 
ciated with strength E and the conductivity will be the 
density dependent factor found before (first used in Eqn. 
I1U|I . Thus, our starting equation is 



dpr 
dt 



d_ 

dx 



J{x) 



d_ 

dx 



D dpr 
dx 



pi 



PrPh 
Ps 



E 



(18) 



As before, we suggest that this driven diffusion equation 
might be obtained from a variation of an appropriate 
large deviation functional, the t = 1 version of which has 
been derived previously [i^, |43 ■ 

As a reminder, our x is actually the fractional length 
along the mRNA, so that it is dimensionless. Similarly, 
the densities are also unitless (e.g., ph € [0,1]), so that 
D/E carries the information of the real length of the 
chain. An estimate of this ratio can be obtained by taking 
the naive continuum limit of a discrete hopping model, 
with the result 1/2N. In the Appendix, we will show how 
this arises from taking such a limit of the MacDonald et 
al. (31 current equation. For here, we may regard D/E 
as a phenomenological parameter. Further, since we are 
focusing only on totally asymmetric processes, we could 
just as well set E to unity (as in Eqn.0), and measure J 
and D "in units of EJ' 

Taking these considerations into account, we seek the 
steady state profile p* {x), as in Eqn. ((T^ . by setting the 
square bracketed terms in Eqn. H18|l to the steady state 
current: 



D{p:r'^-p;pi{p: 



-J* 



In terms of our "natural" variable x (Eqn. 113(1 , this equa- 
tion reduces to 



X 



dx 



-1 



D[l + i£-l) x] 



JX + J*-X{1-X) , (19) 



where J = J* {£ — 1). Thus, we see that this £ > 1 
generalization is quite similar to the £ = I equation: p' = 
[p{l-p)-J*]/D. 

Now, the zeros of x' will play an important role, oc- 
curring at 



I- J + R , 1- J -R 
X> = 7^ and x< = 7\ , 



(20) 



where 



R = \l [1- jY -iJ* 



It is not surprising that the maximal current J (i.e., 
^1 -|- v^j in our system of units) is a key player, so 



12 



that R can be written in the following form: 



J* - J (J* - J) 



Here, J = (l - Vl) appears as a natural counterpart 



to J. Another advantage of this form is its relationship 
to the extremal current principle. As we will see, for 
finite systems, J* will be slightly (O (l/iV^)) larger than 
J for the maximal current phase, giving us complex roots 
and profiles with infiections. On the other hand, for the 
other two phases, J* < J, leading to fixed points in the 
iV — !■ cx) limit. 

Eqn. H19() can be integrated to find an implicit func- 
tion for the density profile in each phase. For explicit 
solutions, we must specify the boundary conditions: 



J* - a[l-p(0)] 





J 



so that 
X(0) = 



! + /?(£- 1) 
p{0)/£ 



piN) 



a - J* 



I- p{o) + p{o)/e a + {e-i)j* 



(21) 



(22) 



and 



x(i) 



PiN) /e 



l-piN)+piN) II 
[l + /3(^- 1)] J* 
- (^ - 1) [1 + /3 (£ - 1)] J* 



(23) 



The boundary conditions can be understood as follows. 
From p- (a) in Eqn. H15(l . we can show that the effec- 
tive reservoir at the initiation boundary has an EPD of 
X_ (a) — a. Eqn. (|14|l then implies the initiation bound- 
ary condition in Eqn. (|21|1 . Particles on the final £ lat- 
tice sites experience no steric hindrance, so the current 
through the final I sites is just 



J* = pr [i) iovi = N 

J* = Ppr (N) . 



1,...,7V~1 



These relations can be used to express p (N) — 

J2i=N-i+i Pr i"^) ill terms of J*, leading to the termi- 
nation boundary condition in Eqn. H21|l . The boundary 
conditions for x '^iH serve to fix the constant of inte- 
gration as well as J*, which is still an unknown in Eqn. 
(I19|l . (The current thus determined is expected to match 
closely the current predicted by Eqn.^l) We discuss the 
various phases separately. 

For smaller values of J* (i.e., J* < J), x{^) has fixed 
points at x> ^'Hd Boundary conditions determine 
whether the steady state profile approaches x< (low den- 
sity) or x>(high density). Each solution corresponds to 
part of the phase diagram shown in Figure |S1 There 
is a "kink" solution when x(0),x(l) G (x<,X>) with 
X (0) ^ x< s-nd X (1) ~ X> that corresponds to the first 
order transition line. 

In the high density phase (/3 < a) where x > X>i the 
density profile is given implicitly by 



[1 + (^ - 1) x>] In -[i + {e-i) X<] In (^^) ^-^ 

I 



(where xo = X (0) is the EPD at initiation) . It may be 
noted that as ^ oo, the bulk density and termination 
boundary density both approach x>- Setting x(l) = 
X> and eliminating J* from the termination boundary 
condition (Eqn. I23|) and the definition of x> (Eqn. I2U|) 
yields a bulk x of (1 - /?) / [1 + /3 (^ - 1)]. We expect 
the bulk density in the high density phase to match the 
termination reservoir density. Indeed, the bulk x value 
calculated here is consistent with the reservoir density 
p+ = 1 — /3 given above (Eqn.^J. 

Though it is impossible to write a closed form for x (x) 
in general, a convenient expression is 



X = X> + (Xo - X>) exp [-p.+x] (24) 

Vxo -x</ 



where 

= J? 

D[l + {£-l)x>] 

and 

^ i-K^-i)x< 
^ i + ie^i)x>- 

Note that i?, x<i and x> can be conveniently estimated 
using the J* value given by Eqn. (|16|l . The advan- 
tage here is that the unknown factor ^ xo-x< ) displays 
only limited variation, since x li^s outside the interval 
[x<,X>]- Thus, we expect to find better and better ap- 
proximations by exploiting an iterative scheme. Starting 
with the first iteration 

X^^^ = X> + (xo - X>) exp [-p.+x] , 
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this procedure involves substituting repeatedly for the x 
in the unknown factor of Eqn. (|24|) . namely, 

r- ' =X>+ Xo - X>) cxp [-tx+x] 

V Xo - x< / 

In Figure^! we show an example of three such iter- 
ations, converging rapidly to the real x{^)i ^iid in the 
inset a comparison with the x {^) observed in simula- 
tions. Near the termination end of the system, the EPD 
from simulations shows a depletion below the bulk den- 
sity. This reduction is characteristic of the high density 
phase and was originally observed in the numerical re- 
sults of 01 . The continuum limit does not capture this 
feature. 

A similar analysis may be carried out for the low den- 
sity phase when x < X<- this case, the boundary layer 
occurs at the termination end, so it is convenient to ap- 
ply the X (1) boundary condition explicitly and use the 
X (0) boundary condition to determine J* . This method 
also works for high and low density phase profiles when 



x(0),x(l) G (x<iX>)j as long as a and fi are not too 
similar. The analysis fails close to the first order transi- 
tion line, as we would expect since Eqn. 1)19(1 predicts a 
"kink" rather than a linear density profile on the transi- 
tion line. 

For the maximal current phase, we know that J* > 
J, so that R is purely imaginary. Since x! is negative 
definite, we have a downward sloping profile, which is 
the generalization of the tan {~x/^) (or cot (x/^)) profiles 
in the ordinary driven lattice gas In the large N 

limit, a macroscopically large region of the profile will be 
almost flat (corresponding to x' ~ 0) , where the density 
assumes the optimal value p (or x — x)- Meanwhile, J* 
will approach the maximal current J from above with 
0{l/N'^) terms. The details are somewhat involved, so 
that we will show the solution and discuss its properties 
only for large N. 

Defining the real quantity R — —iR, Eqn. H19|l can be 
integrated to show that the steady state density profile 
is given by 



In 



x'-[i-J X + J 



+ 2(^- 1) (l - j) 



X'o- i^-J Xo + J* 



Rie-1) 



%) = - 



2x 

D {£-!)' 



(25) 



where 



here can be written as arccot 



6 (x) = arctan 



R 



1-J-2X 



and Xo = X (0) is given by Eqn. 1)22(1 . Evaluating Eqn. 
((25(1 at the termination end {x = l,x = x(l)) ^^d us- 
ing Eqn. 123() . we arrive at an equation for determining 
J* in terms of the control parameters a, fi. Needless 
to say, this equation is too complex to solve analyti- 
cally. Nevertheless, we can gain some insight by con- 
sidering the large N limit. Defining e = 1/N for conve- 
nience, recall that D is 0{e). Next, let us assume that 
J* = J + O (e^) and show that it is justified later. This 

leads us to R being O (e) and (l - /2 = x + O (e^) , 

so that x>,< = X + O (e^) ± iO (e). Thus, the quadratic 
form in the argument of the In in Eqn. I(25() never be- 
comes smaller than R^/4 — O(e^), so that this term 
never exceeds O (In A^). Since the other terms are gener- 
ally O (N) (from D and i?), we can write an approximate 
equation by neglecting the In term, leaving an expression 
for 6, and taking the co-tangent of both sides: 

X - X = cot [x/^ + const.] + O (e^) , (26) 
where ^ = 2DVf./R is of O (1). Note that the constant 



2 (xo — x) / R , which van- 
0. Thus, for typical 



In 



of Eqn. 123) ■ As Figure ITTI shows. 



ishes as R/ (xo — x) in the limit e 
values of x, we see that x = X+O (e). For the termination 
end, J* must be carefully fixed so that the argument of 
the co-tangent in Eqn. ((2(j() approaches n in an appropri- 
ate way for the right hand side to equal x (1) ~ X- This 
requirement is consistent with the original assumption 
that J* = J + 0{e^). 

Returning to the finite A^ case, although x (x) can- 
not be determined explicitly, we find that an iterative 
method is again helpful and converges rapidly. This time, 
successive iterations of x sue substituted into the term 
x"-(i-J)x+J' " 

xg-(l-j)xo + J* 

there is good agreement between the steady state pro- 
files from simulation data and this mean field theory. 



IV. DISORDERED HOPPING RATES IN THE 
OPEN SYSTEM 



We finally turn to the difficult problem of the TASEP 
with extended objects and quenched disorder. At each 
codon of the mRNA, a ribosome translating the mRNA 
must wait for an appropriate aa-tRNA to decode that 
codon. aa-tRNA availability ranges over about an order 
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FIG. 10: Rapid convergence of iterative method to x{^) for small x. (Larger x values are omitted because all iterations of x 
are nearly identical.) x''^'') x'^'i ^-nd x'''' ^^re represented by squares, triangles, and x's, respectively. The system considered 
here is a = 1, /3 = 0.13, A'^ = 1000, £ = 12. D was set to 1/3N to obtain a good fit. Inset compares the actual steady state x 
from simulations (bold curve) with the predicted x'''' (lighter curve). In the profile from simulations, particle depletion near 
the termination end can be seen. 
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FIG. 11: Actual and predicted steady state density profiles 
in the maximal current phase. The system considered here is 
a = 0.5, 13 = 0.5, N = 800, £ = 12. D was set to 4/5N to 
obtain a good fit. J* was set to J + 3.24 x 10~^ to satisfy 
the termination boundary condition. Profiles are the actual 
steady state x from simulations (bold curve) and the pre- 
dicted x'^'' (lighter curve). 



of magnitude 45], and it is thought that the elongation 
(hopping) rates may also range over an order of magni- 
tude, with the slowest ones comparable to the initiation 
and termination rates [4^ . 

Since previous studies of the TASEP with quenched 
disorder involve £ = 1, their relevance to the process 



of translation may be questioned. Indeed, we believe 
that the system will display essential differences when 
these point particles are replaced by extended objects. 
First, as we have seen in the uniform elongation studies 
above, generalizing to the £ > 1 case led to significant 
changes. Now, with non-uniform rates, steric hinderance 
should play an even larger role, since extended particles 
can block multiple sites with potentially different rates. 
Thus, we will devote this section to a limited study of 
the effects of £ > 1 on open systems with non-uniform 
elongation rates: (a) systems with a single slow site, (b) 
bounds on the current in general open systems, and (c) 
as an illustration of systems with full disorder, simulation 
results for a real mRNA sequence. 



A. Simple model of a single internal blockage 

We begin with the simplest form of quenched disorder, 
namely, having a single internal site i with a reduced 
elongation rate: r < 1. If this blockage were moved to 
the termination site, then r would carry the label /3 to 
conform with the notation above. In most of the simula- 
tion studies for this section, we focus on a — (3 — 1 with 
the blockage at the center of the mRNA (i.e., i — 100 
with N = 200). Though the naive expectation is that 
the blockage should limit the current in the same way, 
regardless of its location, we observe that the current 
is noticeably reduced if the slow site is internal. This 
reduction can be understood as follows. Once a particle 
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moves past the slow site, the particle behind it is not nec- 
essarily free to move to the slow site. Instead, because 
of the random sequential updating, the particle behind 
has a nonzero chance of being blocked by the particle 
ahead. With just two particles in the system, this effect 
can be accounted for exactly. Beyond the scope of this 
paper, the details will be published elsewhere. Here, let 
us present some simple heuristic arguments which lead 
us to results that agree well with simulation data. 

Extrapolating from the current in the uniform system 
fEan.[TB|l. we speculate that 



J 



Ps 



free 



(27) 



where Psource is the average coverage density behind the 
slow site and Pfree is the approximate probability for a 
particle just beyond the slow site to advance. Lastly, 
T is the average time to travel through the slow site 
and the i — 1 sites that precede it. For example, in the 
termination-limited case, psource = ^ — fi,T = + (3, 
and Pfree — 1 (no steric hindrance beyond the termina- 
tion site), giving J ^ (3 {1 - (3) / [1 + (3 {£ - 1)]. We esti- 
mate Psource from the bulk density induced if the lattice 
were truncated immediately after the slow site. Pfree 
can be estimated from the density-dependent head spac- 
ing in a closed system (Eqn. ISJ, using the bulk density 
that would be induced if the lattice began with the slow 
site. 

To determine T approximately, we consider the be- 
havior of two particles in an infinite system with a single 
slow site (located at the origin). Denoting the position 
of the leading and following particles by ^ and 77, respec- 
tively, the evolution of this system can be regarded as a 
random walk in the ^-77 plane, confined to ^ — rj > £. Dur- 
ing each time step, the walker moves, with equal prob- 
ability, either "upwards" (?7 — > ?7 + 1) or "to the right" 
(C C, + 1). When it arrives at the ( = r] + £ line, 
the walker remains stationary with probability i (and 
moves "to the right" otherwise). Assuming that there 
are m + 1 holes between the particles immediately after 
the leading one leaves the slow site, so that the walker is 
"initially" located at (Cv) = (1;— ^ — "^)j we computed 
numerically r™, the average time for the walker to arrive 
at the T] = line (i.e., for the second particle to reach 
the slow site). Note that these r's take into account the 
steric hindrance due to the leading particle, which we 
assume is always free to move. To extract T, we make 
the assumption that the left particle advances its first 
m + 1 steps in time m + 1, leaving £ — 1 more steps to 
reach the slow site. Now, the probability p (m) for find- 
ing a gap of TO holes before the leading particle passes the 
slow site can be estimated by inserting the bulk density 
before the slow site (psource) in Eqn. Taking the av- 
erage over these initial starting positions and accounting 
for the time to move over the slow site (l/y), we obtain 
P ~ Sm I"^™ ^ ("t + l)]p (to) -I- 1/r. Predictions for the 
current from Eqn. (|27|l . compared with the actual cur- 
rent from Monte Carlo simulations, are shown in Table 



TABLE I: Actual and predicted reductions in current due to 
an internal slow site with rate r. Currents with the slow site at 
one end are listed for comparison. All currents are multiplied 
by 100. Simulations were performed with A'^ — 200, a = (3 = 1 
with the blockage at the center. These results are statistically 
the same if the blockage is placed at site 20 or 180. 
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Tmax 
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0.01 


13.5 


0.89 


0.86 


0.87 


0.1 


12.9 


4.29 


3.87 


3.90 


0.2 


12.3 


5.00 


4.68 


4.57 



^ for several values of the slow rate r. Despite the ap- 
proximations involved in this approach, the agreement is 
surprisingly good. 



B. Bounds for the current 

To model a real mRNA, we must allow for arbitrary 
translation rates associated with each codon. Let us de- 
note the rate at codon i by ki. Due to the excluded 
volume constraint, however, it is meaningful to consider 
also Ki^i, the maximum rate for a ribosome to translate 
a stretch of £ sites beginning with site i: 



Now, consider a "window" of any stretch of i consecutive 
sites in the lattice. If one particle is moving within this 
window, the following particle must wait until the first 
one passes through entirely before it can begin translat- 
ing these t sites. Thus the characteristic time associated 
with the current must be at least the time required to 
translate the slowest stretch of £ codons. In this way, we 
find an upper bound to the current, i.e.. 



J < min K£ 
ie{i.,...,N-e+i} 



(28) 



One might imagine the slowest segment of £ codons act- 
ing as a "gatekeeper" and preventing the current from 
exceeding the value in Eqn. it^ . 

To arrive at a lower bound, we need only replace the 
elongation rates at each site by the slowest elongation 
rate. From Eqn. Ijltill above, we have the current of a sys- 
tem with uniform rate unity. Thus, the minimum current 
for the disordered system is simply 



J > ( min ki J (a, (3) . 

ie{l,...,Ar-l} ' 



(29) 



where a is the ratio of the initiation rate to the slowest 
elongation rate, and similarly for (3. Though the gap 
between these bounds for a real system may be too large 
to be of significant predictive value, they can provide 
some guide to our understanding of the current. 
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FIG. 12: Steady state coverage density p (upper curve) and 
maximum translation rate 7^12 in each window of ^ = 12 
codons (lower curve) for the ompA gene of E. coU when the 
elongation rates are limiting. Elongation rates at each codon 
were assumed proportional to availabilities of corresponding 
tRNA. 



C. Simulation of a real gene sequence 

To illustrate the full problem of disorder, we have 
simulated translation of several real mRNA sequences 
from Escherichia coli strain MG1655, obtained from [47| . 
Elongation rates at each codon were estimated using 
commonly accepted values for the availability of tRNA 
in E. coli 45] . The rate at each codon was assumed pro- 
portional (with an arbitrary proportionality constant) to 
the availability of the tRNA decoding that codon, as in 
. Corresponding data were not available for estimat- 
ing initiation and termination rates, so a range of initi- 
ation and termination rates was studied. We assumed 
that ribosomes cover i = 12 codons 

Figure^! shows the steady state coverage density pro- 
file for the reasonably well studied gene ompA when the 
elongation rates are limiting. The lower curve in the fig- 
ure shows the maximum translation rate K12 for each 
window oi I =12 codons, calculated from 



K 



12 



IT 

q—l ^ 



for each site i. The minimum of the K12 values is 4.51 
in arbitrary units, thus giving an upper bound for the 
current (Eqn. 128(1 . This value is significantly higher than 
the actual current of 3.52. (For comparison, Ean.l29lgives 
a lower bound of 0.99 for the current in this system.) It 
should be noted that the minimum K12 occurs at codon 
71, which is approximately the location behind which the 
ribosome density is very high, due to ribosomes "piled 
up" behind the slow region. In general, lower values for 
the rate K12 correspond to higher ribosome densities, 
and higher K12 to lower ribosome densities, leading to 
an approximate symmetry between K12 and p. Thus 



the K12 values are useful in understanding the ribosome 
density profiles observed. 



V. CONCLUSIONS 

This work generalizes the well-studied £ = 1 TASEP 
model to particles with extended sizes. Since there is 
a difference between particle density pr and density of 
occupied sites p, the familiar particle-hole symmetry 
{p 1 — p) takes the form pr ^ 1 — p here. Exact 
results for TASEP on a uniform ring, including probabil- 
ity distributions for the current and for particle headway, 
were found. Particularly useful in the latter part of this 
study is the new current-density relationship, which is 
interpreted as a novel density dependent mobility factor 
(or "diffusion constant" ) . 

An extremal principle |3^ based on domain wall theory 
allowed the closed system current-density relation to pre- 
dict currents and bulk densities in the uniform open sys- 
tem as functions of the initiation and termination rates a 
and (3. The phase diagram for the open system was thus 
determined using the extremal principle. Domain wall 
theory also provided an explanation for the linear den- 
sity profiles and other unique characteristics observed at 
the first order transition between the high and low den- 
sity phases. Given the ability of this theory to describe 
particles with length ^ > 1, it might be exploited further 
to determine fluctuations in number of bound particles 
in the steady state and behavior in the pre-steady state 
regime, as has been done for the £ = 1 system 49]. 

Based on the new mobility factor, a simple contin- 
uum limit led to a differential equation for the density 
profile in open systems. Though non-linear in general, 
this equation can be transformed, in the case of a sym- 
metric exclusion process, to the familiar linear diffusion 
equation for an effective particle density x- As a result, 
the stationary density profile, expressed in terms of x, 
is again linear. For the asymmetric exclusion process, 
the differential equation is more complex, due to an ex- 
tra term to account for the drive. Analytic expressions, 
albeit implicit, for the stationary profiles were obtained 
and solved numerically through a rapidly converging, it- 
erative procedure. The results matched simulation data 
closely. However, the predictive power of this approach is 
limited, since the "coarse-grained" parameter D cannot 
be derived from the microscopic jump rates. Naive esti- 
mates for it are of the right order of magnitude but quan- 
titatively inaccurate, reflecting the importance of particle 
correlations. Regarding D as a phenomenological quan- 
tity, we simply fit it to data. At a more detailed level, 
open questions about density profiles remain. In partic- 
ular, our continuum limit cannot explain the particle de- 
pletion and period-^ structure near the termination end 
of high density systems. Profiles of the particle density 
Pr are also of interest, though beyond the scope of this 
simple continuum theory. Absent from the £ = 1 TASEP, 
peaks with spacing £ extend far behind a blockage, with a 



17 



decay length well beyond typical microscopic scales (data 
not shown). Evidently, when extended objects are in- 
cluded, even a simple uniform TASEP tantalizes us with 
a rich variety of novel behavior. 

Finally, systems with quenched disorder in the particle 
hopping rates were briefly considered. Effects of a single 
internal slow site on the current were estimated with fair 
accuracy by considering the average delay a particle near 
the slow site experiences due to a particle ahead of it. 
Also, bounds on the current were determined for general 
disordered systems. A real gene sequence was simulated, 
leading to a complicated density profile. The parameter 
Ki^i, the maximum rate to translate a stretch of i sites 
beginning with site i, proved helpful in understanding the 
shape of the density profile, but the disordered system 
remains far from solved. 

We close with some speculation about the relevance of 
this work to an understanding of translation. The ad- 
vent of functional genomics technologies to measure si- 
multaneously mRNA and protein expression profiles from 
many thousands of genes provides special opportunities 
to begin to understand gene expression regulation. Our 
results for currents (i.e., protein production rates) and ri- 
bosome densities for uniform systems cannot be directly 
compared with such experimental data from typical bac- 
terial cells because translation is not approximated well 
enough by a uniform system. However, it is possible to 
use the uniform system results to interpret data from an 
mRNA artificially constructed to be uniform. Although 
there are no reports in the literature of systems that are 
approximately uniform, the use of an in vitro transla- 
tion system |50| provides an opportunity to make such 
experimental observations. 

It is known that the relationship between mRNA and 
protein levels in typical cells is nonlinear |ll|,il2|- Specifi- 
cally, when cells are grown under two different sets of con- 
ditions, the amount of protein corresponding to a partic- 
ular gene may be down-regulated while its corresponding 
mRNA is up-regulated, and the opposite may be true for 
other genes measured from the same samples. We expect 
that in biological systems, the initiation rate a should 
be an increasing function of the availability of ribosomes 
within the cell. The protein production rate (current 
J, Eqn. I16|l is, in turn, a non-decreasing function of a. 
This analysis thus suggests that the observed nonlinear 
relationship can arise from changes in the availability of 
ribosomes given the nonlinear relationship between J and 
a. However, this situation would cause all protein pro- 
duction rates to change in the same direction, that in 
which the ribosome availability changes. It would not 
permit some proteins to be up-regulated while others are 
down-regulated — which is the situation observed exper- 
imentally. Further nonlinearity would arise if mRNA's 



were to compete for available aa-tRNA as well as for ri- 
bosomes. Detailed modeling of this effect will require a 
better understanding of systems with quenched disorder, 
in which the elongation rates result from aa-tRNA avail- 
ability. 

Acknowledgments 

We thank Johannes Hager, James Sethna, Beate 
Schmittmann, Vassily Hatzimanikatis, Amit Mehra, Tom 
Chou, and Anatoly Kolomeisky for their helpful sug- 
gestions. RKPZ acknowledges support from the Na- 
tional Science Foundation through grants DMR-9727574 
and 0088451. KHL acknowledges support from the NSF 
through grants BES-0120315 and 9874938. LBS was sup- 
ported by an NSF Graduate Research Fellowship and a 
Corning Foundation Fellowship. This research was con- 
ducted using the resources of the Cornell Theory Cen- 
ter, which receives funding from Cornell University, New 
York State, federal agencies, foundations, and corporate 
partners. 



APPENDIX: DERIVATION OF MODIFIED 
DIFFUSION EQUATION FROM DISCRETE 
VERSION 

A modified diffusion equation qualitatively equivalent 
to Eqn. (|19|l can be obtained by taking the naive con- 
tinuum limit of the discrete mean field equations of Mac- 
Donald et al. We begin with Eqn. (7) of 7] for the current 
from lattice site j to j + 1: 



(0) 



and make the following correspondences with our contin- 
uum notation: 



■f 

(i) 

3 

L 

q 



Ph {x) 
Pt [x) 

I 
J 



where x lies in [0, N] here. Thus we have 

Pr {x) Ph {x + 1) 



J{j^J + l) = 



Ph{x+l)+pr ix+£) 



(A.l) 



Performing a series expansion of Eqn. (jA.l|) and keeping 
terms up to second order, we find that 
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by using ph {x) = 1 — ipr {x). Similarly, 
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It can then be shown that 



equations of predict 



d T 

— = J(j^j + 1)- J(j-l^j) 



d_ 

dx 
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( PrPh \ 
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\ Ps J 



This leads to the modified diffusion equation 
'1 + £{!-£) pI 



J = D 



dpr 

dx 



E^lEh. , (A.2) 

Ps 



where we have defined D — 1/2 and E — 1. 

Finally, we rewrite Eqn. 1|A.2|I in terms of the effective 
particle density x- Thus we find that the mean field 



dx 
dx 



l + {£-l)x 



i + ii~i)x-{e-i)x' 



(A.3) 

We expect Eqn. (|A.3|I to be comparable to Eqn. 
the steady state density profile equation derived previ- 
ously from simple arguments. Indeed, both equations 
give the same fixed points for x, thus producing quali- 
tatively identical density profiles. Further, their quan- 
titative differences have little effect on the shape of the 
density profile (data not shown). 
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